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The powerful molecular dynamics (MD) simulation is basically based on a picture that the atoms experience 
classical-like trajectories under the exertion of classical force field determined by the quantum mechanically 
solved electronic state. In this work we propose a quantum trajectory approach to the MD simulation with 
surface hopping, from an insight that an effective "observation" is actually implied in the MD simulation through 
tracking the forces experienced, just like checking the meter's result in the quantum measurement process. This 
treatment can build the nonadiabatic surface hopping on a dynamical foundation, instead of the usual artificial 
and conceptually inconsistent hopping algorithms. The effects and advantages of the proposed scheme are 
preliminarily illustrated by a two-surface model system. 



A full quantum mechanical treatment for molecular dy- 
namics (MD) would break down along the increase of 
atomic degrees of freedom. As a result, in practice the 
MD technique based on an assumption that atomic mo- 
tions are governed by classical mechanics has proved to 
be a very powerful tool [0-^]] . The MD simulations usu- 
ally rely on the Born-Oppenheimer (BO) approximation, 
which decouples the electronic and atomic motions. That 
is, the atoms experience classical-like trajectories under 
the exertion of appropriate classical force field, which is 
determined by a single potential energy surface (PES) as- 
sociated with a single electronic state. On the other aspect, 
the electronic states are solved from the time-independent 
Schrodinger equation for a series of atomic geometries 
(configurations). 

In many situations, however, when the energy separa- 
tion of different PESs becomes comparable with the mag- 
nitude of the nonadiabatic coupling (typically in the prox- 
imity of conical intersection, see Fig. 1), the BO approxi- 
mation may often fail. In this case, since each nondegen- 
erate electronic state defines a distinct BO PES, a transi- 
tion between electronic states would change, often drasti- 
cally, the forces experienced by the atoms. Proper account 
for this nonadiabatic effect is of crucial importance in 
practice. The treatment of nonadiabatic effects in MD has 
a long history. The most widely applied include the so- 
called "Ehrenfest" or "time-dependent-Hartree mean-field 
(TDHMF)" approach [ffl Pj-fl K the "trajectory surface- 
hopping (TSH) methods [pH16||, and their mixed schemes 
[[i7|-|r^]. The former approach views that the electronic 
wavefunction might in general be a linear combination 
of the BO adiabatic functions, the atomic effective poten- 
tial is thus calculated by averaging the electronic Hamil- 
tonian over such wavefunction, while the BO adiabatic 
functions are determined self-consistently with the atomic 
trajectory. The TSH, a more appropriate methodology for 
dynamics propagation of nonadiabatic systems, is how- 
ever based on an insight that the trajectories should split 
into branches, i.e., each trajectory should be on one state 



or another, not somewhere in between. And, the trajec- 
tories distribution is achieved by allowing hops between 
surfaces according to some probability distribution. 



(A) 




(B) 



Reaction Coordinate (R) 



KV\AM 





'Electronic address: ixinqi@bnu.edu.cr 



Electronic States Atomic Potentials 

FIG. 1: (A) Schematic atomic potentials in terms of adiabatic 
(solid) and diabatic (dashed) representations for the electronic 
ground and first excited states along the reaction coordinate. 
(B) Measurement analogy by dividing the molecular system into 
two subsystems, and regarding the atomic part as a "measure- 
ment apparatus" which continuously probes the distinct elec- 
tronic states by the experienced distinct forces. 

A variety of trajectory-based methods have been devel- 
oped. Among them, the most typical one is Tully's/ewesf 
switches surface hopping (FSSH) algorithm [[l5|]. This al- 
gorithm is based on the following insights/arguments: ( i) 
The atomic trajectories determine the probabilities of elec- 
tronic transitions and the electronic transitions, in turn, 
strongly influence the forces governing the trajectories. 
(ii) The atoms evolve on individual single PES and the 
nonadiabatic effects are included by allowing hopping 
from one PES to another according to the weight of the 
respective electronic state. In the past two decades, stim- 
ulated by Tully's work, a variety of variants of the FSSH 
algorithm and a large number of applications have been 
carried out [p 



2 



Nevertheless, the coherence between states, maintained 
in the FSSH algorithm is usually wrong because of the in- 
dependent trajectory approach (see the recent review [||]). 
The underlying conceptual difficulty can be further un- 
derstood as follows. In the FSSH algorithm, the nonadia- 
batic effect is treated by allowing hopping from one PES 
to another, with the hopping probability determined by the 
weight change of the respective electronic states in quan- 
tum superposition. Physically, however, the atomic clas- 
sical trajectory along a specific individual PES implies 
that it must collapse the electronic state onto the corre- 
sponding BO wavefunction, according to the distinct clas- 
sical force experienced. In other words, we can no longer 
treat the electronic state in a quantum superposition of the 
many different components, as in contrast proposed in the 
FSSH approach [Oj. As shown in Fig. 1(B), through an 
analogy with quantum measurement or with the more pop- 
ular Schrodinger's Cat paradox, the FSSH algorithm sim- 
ply means that we are still treating the radioactive decay 
in a quantum superposition state, while having found the 
Cat definitely "alive" or "dead". 

In this work, rather than designing certain wiser algo- 
rithms for the TSH probability, we would like to view the 
problem in an alternative way and accordingly propose a 
quantum trajectory approach to the MD simulation. That 
is, we are able to develop a self-consistent in physics and 
hopefully highly implementable stochastic MD scheme, 
which quite naturally renders the TSH in the most funda- 
mental notion of quantum jump. 

Formulation. — For the purpose to be clear soon, we 
introduce the electronic Hamiltonian by subtracting the 
atomic kinetic energy from the total Hamiltonian as 
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Here, the first term describes the kinetic energy of elec- 
trons. We use r denoting all the electronic coordinates, 
i.e., r = {rj\ j = 1, 2, • • ■ }. And, formally, R denotes the 
atomic configuration. 

We expand the electronic wavefunction by a set of 
known basis functions, ^(r, R) = ^ - Cj(t)(j>j(r, R). In 
principle, {<f>j(r, R)} can be any electronic basis func- 
tions, while in practice it would be convenient to chose 
them as the Born-Oppenheimer adiabatic wavefunctions. 
That is, H e i(r,R)(f>j(r,R) = £j(R)0j(r,R), which 
means that (f>j (r, R) are the instantaneous eigenstates of 
H e i(r,R), for a given R = R(t). These basis wave- 
functions, together with the corresponding eigen-energies 
Ej (R), are the standard output of the usual quantum chem- 
istry computation. Using the selected basis functions, we 
further introduce 

V jk (R) = (hirMlH^fa^RWMr,*)), (2) 
d jk (R) = (^(r.RJlVR^fcCr.R)). (3) 

Obviously, for the BO adiabatic basis, Vj k (R) is diagonal, 
V jk (R) = ej(R)6 jk = Vj(R)5 jk . Each Vj(R) is known 



as the BO potential energy surface (PES). The quantity 
djfc(R) characterizes nonadiabatic coupling between the 
( jth and fcth) PESs, which has an effect of driving surface 
hopping. This can be seen more clearly by the follow- 
ing. The coherent electronic evolution is governed by the 
Schrodinger equation itffil = H e i^, which results in a set 
of coupled equations for Cj (t) 

ih6 i = Yl l v M R ) ~ iff* ' <M R )] c fc 

k 

= Y^H ]k (R)c k . (4) 

j k 

In this expression, we can understand Hj k (R) as an ef- 
fective Hamiltonian matrix in the selected basis, which 
straightforwardly describes the evolution of the state 
amplitudes Cj(t). In terms of a projective operator 
form, the effective Hamiltonian can be expressed as 
H(R) = J2 jk H jk m^(R)){MR)\- Using the prop- 
erty d* fc (R) = — dfcy(R), one can prove that H(R) is 
Hermitian, as it should be. Moreover, corresponding to 
the BO adiabatic basis, we have even simpler results: 
HjjiR) = £j (R); and H jk (R) = Q jk (R) = -ihR ■ 
d 3 fc(R), if j 7^ k. In this way, we see that it is the prod- 
uct of R and dj k (R) that characterizes the nonadiabatic 
coupling between the BO potential surfaces. 

It was merely based on Eq. (|J), which describes the 
quantum mechanical evolution of the superposition ampli- 
tudes of electronic states, that the highly celebrated FSSH 
algorithm was proposed [15|, with the following central 
idea. From the "normalized" population change of the 
electronic state \(j>j(R)), defined as pj(t) = [\cj(t)\ 2 — 
\cj(t + At)\ 2 ]/\cj(t)\ 2 , one performs a standard Monte- 
Carlo choice (with this probability Pj(t)) to determine 
whether or not a hopping event should take place, from 
the BO potential surface Vj (R) to another one. 

Intuitively, this looks indeed a promising algorithm, 
since the BO potential surface t^-(R) is anyhow associ- 
ated with the electronic state \<fij(R)), one can then con- 
clude that the change of the electronic occupation proba- 
bility must imply a surface-hopping event occurred. Nev- 
ertheless, in each single realization of MD trajectory, the 
atomic motion experiences a series of distinct PESs (but 
not their weighted superposition), together with succes- 
sive hopping between them. Based on the entanglement- 
type correlation between the electronic states and the 
atomic PESs, see Fig. 1(B), the classical trajectory treat- 
ment for the atomic motion must collapse the electronic 
state to an individual one, with the result determined by 
the classical force experienced by the atoms. In other 
words, we can no longer inappropriately maintain the 
electronic state in a quantum superposition of two (or 
many) different components \<fijJR)), as in contrast as- 
sumed in the FSSH approach [n5p. In essence, using the 
term of Schrodinger's Cat paradox, the mixed quantum- 
classical FSSH treatment is equivalent to treating the ra- 
dioactive decay still in "quantum superposition", while 
one has definitely found the Cat "alive" or "dead". 



3 



In order to eliminate the above inconsistency in the 
FSSH algorithm, we propose that one should apply a 
measurement-based description to the electronic state evo- 
lution, instead of using the R-dependent Schrodinger 
equation. Indeed, since R enters the Schrodinger equa- 
tion, this partly accounts for the effect of the atomic mo- 
tion on the electronic state evolution. But this is not 
enough. Actually, the atomic motion is continuously get- 
ting the state information of the electrons, via the corre- 
lation between the PES and the electronic state. As a re- 
sult, in addition to involving "R" as an external parameter 
in the electronic Schrodinger equation, one should at the 
same time account for the backaction of the electronic- 
state information gain. This means that the atomic mo- 
tion is in essence making a continuous quantum measure- 
ment on the electronic state. The distinct "force" experi- 
enced by the atomic trajectory motion, effectively, plays 
the same role as the meter's output in quantum measure- 
ment process. 

After the above concept ual preparation, we can then ap- 
ply the established quantum trajectory equation (QTE) to 
account for the backaction effect owing to state informa- 
tion gain as follows 

p = ~[H(R),p]+^r w 2?[M 3 -(R)]p 

3 

■ Y^\~n MAR- p^(t). (5) 
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Here, p is the electronic state density matrix, with ele- 
ments pjk- The diagonal elements describe population 
probabilities on the BO states, while the off-diagonal el- 
ements describe quantum coherence between them. Note 
that, on the left hand side of Eq. (0), we have made a con- 
vention that p represents J2jk Pjkl&j W) W l> ^ ut not 
iiEjkPM^KMm- Accordingly, H(R) in the 
above QTE is the effective Hamiltonian defined in Eq. 

In Eq. (Q), the measurement operator, Mj (R) = 
\<pj(R))(<f>j(R)\, corresponds to an identification (dis- 
crimination) to the electronic state <fij(R) by the atomic 
motion. It is well known that any state discrimination 
will destroy the quantum coherence (the quantum super- 
position). This effect is described by the second term in 
Eq. (Q), with the dephasing rate corresponding to dis- 
crimination of the state 4>j (R). Here the Lindblad-type su- 
peroperator means, V[Mj]p = M,p.\l ] - \{M]Mj,p\. 
Finally and very importantly, the last term in Eq. (ph 
accounts for the effect of state information gain. The 
involved superoperator, more explicitly, is defined as 
H[Mj}p = M jP + pM] - (Mj + M])p, where (Mj + 

M]) = Tr[(Mj + Mhp\. £j(t) are the Gaussian 
stochastic noises, satisfying the ensemble average prop- 
erty E[£j(i)ffc(i')] = S jk S(t-t'). Notice that the last 
term in Eq. ^ is not originated from any external noise, 
but from an intrinsic stochastic collapse (quantum jump) 
associated with "observation", i.e., the quantum measure- 
ment postulate. 

Applying Eq. (^) in a strong observation regime (with 



large T^f), one can generate the desired TSH picture. That 
is, the large T^j correspond to MD simulation under dis- 
tinct forces. Then, under the influence of the atomic tra- 
jectory motion, the nonadiabatic coupling will drive a se- 
ries of stochastic jumps between the BO states. This is 
nothing but the desired "surface hopping" behavior which 
can be generated, in the existing schemes, only by con- 
structing various hopping "algorithms", non-dynamically. 

Alternatively, applying Eq. (|j) in a weaker observa- 
tion regime (with smaller Y^j), while the atomic mo- 
tion propagates along a single PES away from the level- 
crossing area, the Ehrenfest-type mean field approach 
is restored in the proximity of the level-crossing area. 
That is, by computing the potential energy via V(R) = 
Tr[H e i(r, R)p(t)}, the effect of quantum superposition of 
the BO potential surfaces is taken into account. In our 
opinion, this effect should exist there. Using the term of 
quantum measurement, this simply corresponds to an in- 
complete observation in the level crossing area where the 
atomic motion, particularly in the case of strong nonadi- 
abatic coupling, does not fully distinguish which poten- 
tial surface, but only experiences a mixed force. We be- 
lieve that applying the proposed approach to the above 
regime is a valuable extension of the hybrid TSH-TDHMF 
scheme. 

Illustrative Demonstration. — Below we illustrate how 
the approach proposed above can generate the TSH be- 
havior, dynamically (or physically in essence), instead of 
using the probability algorithms to generate the "hopping" 
events. For simplicity, we restrict our demonstration to 
a one-dimensional (ID) two-surface example. Because 
of the ID nature, in what follows we denote "R" sim- 
ply by "i?". To be specific, we model the two BO po- 
tential surfaces as schematically shown in Fig. 2(A). For 
\R — i?o| > a, we assume E2(R) — £i(R) = &\R — R \; 
while for \R — Rq\ < a, the adiabatic energy difference 
is a constant, £2{R) — £i{R) = A. Here, the scale of the 
level-crossing area is determined by a = A/a. Rather 
than a real MD simulation, following the Landau-Zener 
model, we assume a constant speed for the atomic mo- 
tion, i.e., R(t) = vt. Moreover, the nonadiabatic coupling 
is parameterized as H\2 = (H21)* = —ihvf3 = —ift£l, 
from an assumption of R ■ di2(R) = v/3. 

In Fig. 2(B) and (C) we display some representative 
quantum trajectories. First, we notice that the nonadia- 
batic coupling (with constant magnitude "f2") can cause 
efficient transitions largely in the proximity of the conical 
intersection area, because of the relatively small energy 
separation between the PESs. Also, in the whole range, 
Eq. (|j) will give a quantum mechanically pure (but not 
unitary) evolution for the electronic state. Fig. 2(B) cor- 
responds to a situation with relatively weak coupling and 
weak observation. In this case, in the intersection area the 
electronic state consists of a superposition of the BO ba- 
sis states. Using this state to calculate the atomic forces, 
desirably, should be an extension of the usual mean-field 
approach. However, if we increase further, as shown 
in Fig. 2(C), in the proximity of the central level-crossing 
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FIG. 2: Demonstration of the surface-hopping behavior. (A): 
The model we used in our numerical simulation. (B) and (C): 
Four representative trajectories (red curves) from two groups of 
parameters, while the cyan curves stand for the results of Vj, = 
in the absence of "observation". Other parameters: A = 1.0 and 
a — 10. In the simulation (also in Fig. 3 and in the main-text 
discussion) we assume T^i = r^2 = IV 



area the sharp "surface hopping" behavior will appear 
eventually. 

We emphasize in this context that, very importantly, 
even for a weak strength the trajectory will gradu- 
ally and stochastically collapse onto one of the energy sur- 
faces. This is in sharp contrast with the prediction of the 
Schrodinger equation. In that case, the state will remain in 
quantum superposition, after the evolution passes through 
the central intersection area. We demonstrate this feature 
by the cyan curves in both Fig. 2(B) and (C), where we 
see that the final state is indeed in superposition. How- 
ever, a continuous "observation" will collapse the state, 
gradually, onto one of the basis states, either \(f>i) or \cf>2)- 
Then the desired picture of trajectory propagation along a 
single energy surface is generated, by this dynamical and 
physical means. 

In practice, to compare the stochastic MD simula- 
tion with experiments, one should ensemble-average the 
stochastic trajectories. For the assumed constant speed 
of atomic motion, the trajectory distribution can be sim- 
ply predicted by averaging the QTE first. This can be 
done easily by removing the last (unraveling) term in 
Eq. (Q) and remaining only the first two terms. We ac- 
tually obtain the usual master equation. Then, from the 
density matrix, we get the probabilities pi(t) = pn{t) 
and P2(t) = P22{t), which give the distribution ratio of 
the atomic trajectories in the asymptotic limit (t — > oo). 
We should note that, however, in practical MD simulations 



one cannot average first Eq. (q|), since for different trajec- 
tories the stochastic atomic forces experienced would be 
different. Obviously, this feature cannot be captured by 
the averaged master equation. In other words, the aver- 
aged master equation does not allow us to correctly com- 
pute the atomic forces, we are then unable to propagate the 
atomic motion based on Newton's law. This is the reason 
that in literature, when one attempts to introduce decoher- 
ence, the master equation approach does not work. 
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FIG. 3: Nonadiabatic transition probability from the initial state 
\<pi(R- 00 )) to the final one \4>2(Roo))- R^oo represent the ini- 
tial and final asymptotic reaction coordinates. This probability 
indicates the distribution of the final "products". (A) Depen- 
dence on the nonadiabatic coupling strength in the absence of 
dephasing (T^ = 0). (B) Dephasing effect on the transition 
probability. Parameters: A = 1.0 and a = 10. 

In Fig. 3(A), we show the dependence of the nonadia- 
batic transition on the coupling strength. We take T^ = 
in order to reveal this mere dependence more clearly. We 
observe a turnover behavior, with an optimal nonadia- 
batic coupling to make the transition reach maximum. 
This feature indicates, in some counterintuitive manner, 
that a faster atomic motion may not necessarily enhance 
the transition probability. This behavior differs somehow 
from the prediction of the Landau-Zener formula, since 
that formula will give a larger transition probability for a 
faster speed. The reason is that, here, in the adiabatic BO 
basis, the faster atomic motion will enhance the nonadia- 
batic coupling between the energy surfaces, while in the 
Landau-Zener model, expressed in a diabatic basis, the 
passage speed across the intersection area has no such ef- 
fect . 

In Fig. 3(B), we show the decoherence effect on the 
nonadiabatic transition. Here we observe a decoherence- 
enhanced transition behavior. But, this feature is not uni- 
versal. It only corresponds to a weak dynamic tunneling 
regime. Actually, not shown in Fig. 3(B), if we tune the 
coupling strength into an intermediate tunneling regime, 
a turnover behavior can appear, whereas a decoherence- 
suppressed transition will take place in strong tunneling 
regime. We also observe that, if an efficient coupling per- 



5 



sists in the intersection area for longer time, the decoher- 
ence will result in a final equal occupation on the both 
energy surfaces, as shown by the curve with il = 1.5 in 
Fig. 3(B). 

Finally, we remark that a decoherence in the TSH 
should arise because the MD simulation is propagated 
along a single trajectory determined by the gradients for 
an specific electronic state. In this propagation, the am- 
plitudes of all other states are artificially restricted to be 
also propagated along the same trajectory, even though the 
gradients of these other states would dissipate these am- 
plitudes along other directions. This effect has been dis- 
cussed in detail and corrections to it have been proposed 
|]3^-^|. For instance, in Ref. [55h, it was shown that the 
divergence between occupation and average population is 
caused by the missing decoherence, and that this discrep- 
ancy can be eliminated if the decoherence corrections are 
incorporated. 
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Appendix A: Measurement Interpretation 

In this appendix we present a further discussion to re- 
late the stochastic MD simulation with quantum measure- 
ment. As mentioned in the main text, we can divide the 
whole system into two coupled electronic and atomic sub- 
systems, as shown in Fig. 1(B). Under the BO approxima- 
tion, the relatively heavy and slowly moving atomic sub- 
system largely experiences a classical-like trajectory mo- 
tion. Therefore, the atomic subsystem looks very like a 
measurement meter, which is continuously "measuring" 
the electronic states by its experienced distinct forces. In 
this context, the PES (or force) is the meter's output (mea- 
surement result), from which the electronic state is in- 



ferred. 

In a more heuristic manner, we may relate the problem 
to the paradox of Schrodinger's Cat. As illustrated in Fig. 
1(B), we assume two electronic states, \g) and |e). Ac- 
cordingly, we have two PESs, V g (R) and V e (R), which 
correspond to the Cat states "alive" and "dead". If we in- 
sist on regarding the whole system as a closed one, we will 
arrive to a (quantum) superposition of V g (R) and V e (R), 
exactly as the superposition of "alive" and "dead" states 
in the Cat paradox. But, as is well known, the Cat para- 
dox can be resolved only by a neglect of the large number 
of microscopic degrees of freedom of the Cat and the sur- 
rounding environment. It is merely this treatment that can 
result in the emergence of classicality which shows the 
Cat in a state of either "alive" or "dead", but not in be- 
tween. We then arrive to two statements: f i) In the semi- 
classical stochastic MD simulation, the essential picture 
of classical trajectory and surface hopping has implied an 
incorporation of certain external environment and an ef- 
fective observation. For a closed quantum system, there is 
no way to generate the jump (hopping) behavior. ( ii) We 
hypothesize that the TSH-based MD simulations are more 
appropriate for many-atom system in condensed phase. 
For an "isolated (clean)" few-atom system, a full quantum 
mechanical simulation should be required. 

In the TSH-MD simulations, we can imagine that the 
"observation" is realized through tracking the forces expe- 
rienced, just like checking the meter's result in the quan- 
tum measurement process. Together with the quantum 
nonadiabatic coupling, successive jumps then appear, ex- 
actly like in the continuous measurement of a driving 
system^lp. Moreover, this approach is quite versatile. 
For instance, it can describe the possible "partial" hopping 
behavior. That is, when the atomic motion passes through 
the conical intersection region, the atoms may not be able 
to resolve the experienced force from which PES. In this 
case, the resultant state is still in quantum superposition. 
However, the state evolution should obey the quantum tra- 
jectory equation, Eq. (Q), but not the Schrodinger equa- 
tion. We believe that this treatment is a valuable extension 
of the hybrid TSH-TDHMF approach [|f7]-(l|]. 
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